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We study a chain oi N + 1 phase oscillators with asymmetric but uniform coupling. This type of 
chain possesses 2^ ways to synchronize in so-called travelling wave states, i.e. states where the phases 
of the single oscillators are in relative equilibrium. We show that the number of unstable dimensions 

of a travelling wave equals the number of oscillators with relative phase close to n. This implies 
that only the relative equilibrium corresponding to approximate in-phase synchronization is locally 
stable. Despite the presence of a Lyapunov-type functional periodic or chaotic phase slipping occurs. 
For chains of length 3 and 4 we locate the region in parameter space where rotations (corresponding 
to phase slipping) are present. 

PACS numbers: 05.45.Xt,87.19.1m,05.45.-a 
Keywords: synchronization, phase slipping 



I. INTRODUCTION 

Investigations of synchronized behavior of coupled 
nonlinear oscillators permeate physics [1, 2], chemistry 
[3, 4], biology [5] and engineering [6]. Dating back to 
Huyghcns, this problem has been explored both mathe- 
matically and experimentally. The groundbreaking work 
of Kuramoto [7] led to the realization that sync;lironiza- 
tion is a ubiquitous behavior in nature and can be ex- 
plained by simple models for interaction between compo- 
nents of a system (see for example the excellent review 
of Strogatz [8]). 

Cohen et al. [9] studied a chain of Kuramoto oscillators 
to explain the so-called fictive swimming observed in the 
Central Pattern Generator (CPG) of the primitive verte- 
brate lamprey. Their purpose was to explain the uniform 
phase lag along the segmental oscillators. Later, Kopell 
and Ermentrout [10] proposed a more realistic model for 
fictive swimming and Williams et al. [11] contrasted this 
model with experimental observations. Studies of the 
CPG led to bio-inspired applications in robotics, for ex- 
ample, the autonomous mobile robotic worm of Conradt 
and Varshavskaya [12] and a turtle-like underwater vehi- 
cle by Seo et al. [13]. 

Cohen et al. [9] observed that in a chain of oscil- 
lators the oscillators at the end points play as special 
role: adapting the natural frequency of the end points 
controls the phase shift between neighboring oscillators 
throughout the entire chain in the phase-locked equilib- 
rium. (All equilibria of the one-dimensional homoge- 
neous Kuramoto model with periodic boundary condi- 
tions were found in [14].) 

In this paper we vary this detuning 5 of the end points 
and the coupling strength (more precisely the ratio k be- 
tween coupling strength down the chain and up the chain) 
to study transitions between phase locked solutions and 
phase slipping. If one introduces the phase differences 
between the oscillators as the new dependent variables 



then the phase-locked solutions are equilibria, and for a 
chain of length A^^-l- 1 there exist 2^ of these equilibria. A 
transition from one phase locked solution to another then 
corresponds to a motion along a heteroclinic connection 
between the corresponding saddles in the phase space. 
The 2^ equilibria can be classified by an index quantity 
Vtt which counts how many of the phase differences are 
equal to tt — ^. This index turns out to be identical to 
the number of unstable dimensions of the saddle (pro- 
vided k > —1) such that connections between equilibria 
of decreasing index i^^ are generic. Solutions with contin- 
uously slipping phases between oscillators show up as ro- 
tating waves in the phase space. Periodic phase slipping 
and its bifurcations can be computed directly if one as- 
sumes, for example, pure phase coupling (coupling of the 
type Ksm{Oj — Oj-i) for neighboring oscillators 6j and 
6*^+1). One of the co-dimension 1 boimdaries of rotating 
waves are non-generic connections between equilibria of 
identical index . We present numerical evidence of this 
for the cases of iV = 2 and = 3 (that is, the case of 
3 and 4 oscillators, respectively). A noticeable difference 
between these cases is that the basin of attraction for ro- 
tating waves appears to be a slightly smaller fraction of 
the phase space for larger N, and the parameter region in 
the (fc, (5)-plane permitting rotations (that is, continuous 
phase slipping) is smaller for the larger N. We conjecture 
that making an oscillator chain longer does not make its 
propensity for phase slipping larger even if the coupling 
strength between neighboring oscillators is not increased. 



II. MODEL DESCRIPTION 

We consider a chain of N + 1 phase oscillators with 

nearest neighbor coupling as shown in Fig. l.This model, 
first proposed by Cohen et al. [9] to explain the 'fictive 
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FIG. 1: Schematic representation of chain coupled oscillators 
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swimming' observed in Lamprey spinal cord, is given by 
d 

for j = 1,2, ...,7V- 1, 

d 



dt 



where 9j and ujj ^ are the phase and natural frequency 
of the j-th oscillator, respectively. The coupling gains 
are denoted by Ky^ and Kj^. The prototypical example 
for the coupling function is T{x) = sinx. This choice of 
the coupling function has the following essential features: 

1. r is 27r-periodic, |r(a;)| < 1, 

2. r(x) has odd symmetry about a; = (that is, 
T{x) = —T{—x)) and even symmetry about x = 
tt/2 (that is, r(7r/2 - x) = r(7r/2 + x)), 

3. r{x) > for X e [0,7r/2]. 



III. PHASE LOCKING AND TRAVELING 
WAVES 

Since the frequencies ojj are non-zero, system (1) does 
not in general have fixed points. If all oscillators are 
identical (that is, Wq = . . . = ujn) Equation (1) admits a 
fully synchronised state 6q = . . . = 0n. Here we focus on 
the so-called phase locked solutions, that is, solutions of 
the form 9i — 9j = const. Of particular importance are 
phase locked solutions where the phase differences on the 
chain are identical, called uniform traveling waves. 

To analyze phase locked solutions we introduce the 
phase differences 



9, {j^l...N), 



(2) 



as new variables, as well as the rescaled frequency dif- 
ferences flj, the coupling ratio k and the rescaled time 

^ncw ■ 

System (1) leads to a system of TV equations describing 
the dynamics of the variables Xj wrt the new time t^cw 

ii = 17i-(l + fc)r(xi)+r(a;2), 
Xi = + kT{xi^i) + T{xi+i) - (1 + k)T{xi), 

i = 2,3,...,7V-1, 
Xn + k^ixN^i) - (1 + k)T{xN), 



(3) 



where Xi (t) is on the unit circle S"'^ . The vector form reads 

x = o + cr(x), (4) 

where x = [xi,X2,-- - ,xn]'^, ft = [fli,fl2,--- ,^nV, 
r(x) ~ [r(a;i), . . . ,r(xAr)], and C is an x iV matrix of 
the form 



C 



-(1 + fc) 1 
k -{l + k) 














-(l + k) 1 
k -(1 + fc) 



(5) 



The fixed points of Equation (4) are given by r(x) = 
—C~^fl (C is invertible for fc 7^ — 1, see Proposition 2) 
and correspond to the phase locked solutions satisfying 
9i — 9j — const of the original system (1). 

If all components of —C~^n are less than 1 in abso- 
lute value, system (4) has 2^ equilibria, because of the 
even symmetry about 7r/2 (r(7r/2 — x) = r(7r/2 + x), see 
Section 2) of the coupling function. A uniform traveling 
wave solution (phase locked solution with identical phase 
differences, i.e. x = [S, - ■ ■ , <5]"^)) exists only if the follow- 
ing conditions on the rescaled frequency differences Q are 
satisfied: 



ni = kr{S), 

a =0, i = 2, 



,N-1, 



(6) 



This means that the frequencies ojj of the original system 
(1) must be of the form 



cjo = + KdT{S), 

Uj = cj, i = 1, ■ 



,,7V-1, 



(7) 



where w e M. In other words, all oscillators must have 
identical natural frequencies except for the two oscillators 
at the boundary. The "detunings" (difference from the 
uniform frequency uj) of the first and last oscillators are 
related to one another via the coupling strengths and 
Kd- 

The two primary parameters affecting the dynamics 
are the coupling strength ratio k and S. Without loss 
of generality we can restrict our considerations to the 
parameter set 

• (5 e [0,7r/2]: only r((5) enters the equation and for 
negative S we can apply the transformation x 1— > —x 
since F is odd. 

• fc > — 1: for fc < — 1 we can apply the transforma- 
tion [xi, . . . , a;jv]ncw = i^Nj ■ • ■ J 2^i]oldi ^now — fcold^j 

and fcnow = — 1/fcoid such that fcnow > 0. Note that 
this transformation involves reversal of time direc- 
tion (since fcoid is negative) so statements on sta- 
bility will have to be replaced by the corresponding 
statements on instability and vica versa. 
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Wc observe that the couphng strength does not en- 
ter Equation (3), which determines the dynamics, at all. 
Only the ratio between down-chain and up-chain cou- 
pling strength matters. The absolute value of the cou- 
pling strength then determines the time scale with re- 
spect to the the original time (Equation (3) is with re- 
spect to a rescaled time Kut). 



IV. STABILITY OF TRAVELING WAVES 

Naturally, wc are interested in characterizing the local 
stability of traveling waves, that is, of equilibria of (4). 
We note that even though the uniform traveling waves 
can occur only in a slightly degenerate parameter set- 
ting, these waves are robust in the sense that slightly 
non-uniform traveling waves will exist for slight pertur- 
bations of these parameters values. For example, non- 
uniform traveling waves exist for a linear gradient fre- 
quency distribution (il = const.). Cohen et al. [9] de- 
rived a necessary condition for the existence of the phase- 
locked solutions. Ermentrout and Kopell [15] showed the 
existence of frequency plateaus when this necessary con- 
dition is violated. 

The linearization of (4) about an equilibrium a;* is J = 
CT'{x^) where C is the coupling matrix given in (5), and 



r'(a;*) =/>-diag(<Ti,...,(T„), 



(8) 



is a diagonal matrix with p = T'{5) > Q (since 5 e 

[0,7r/2]) and Ui = ±1, depending on whether the ith 
component of a;* is equal to J or tt — (5. Even though the 
eigenvalues of the matrix C are known, i.e. 



Aj = -(l-hA;)-|-2\/fccos 



N + V 



1, 



(9) 



the eigenvalues of the product Cr'(a;*) are not known 
analytically. Therefore, the goal of this Section is to es- 
tablish a simple connection between the number of stable 
and unstable directions of an equilibrium a;* (and thus 
its stability) of (4) and the number of components of a;* 
equal to 5 and -k — 5. To establish the result, we associate 
the four indices with each equilibrium 

vq{x^) = number of components of equal to 5, 
2^77(2;*) = number of components of a;* equal to tt — (5, 
^u{x*) = dimension of unstable manifold of a;*, 
i^six*) = dimension of stable manifold of a;*. 

(10) 

The relation between these indices is given by the follow- 
ing 

Theorem 1 (Dimension of invariant subspaces). Let x^, 
be an equilibrium of (4). If k > —1 then 

fo(a;*) = i^six*) and ^',r(a;*) = i^u{x*)- (11) 



Proof Relation (11) is proven indirectly via the 
Lyapunov-type functional 



^ r 

E{x)=Y, / r{y)y-T{S)xi 
i=i L^o 



(12) 



We observe that the time derivative of E along trajecto- 
ries of (4) is 



N 



^ (^(r(xi) - mf + (r(x^) - r(<5))^ 

JV-1 \ 



+ J2{T{x,)-rix.+,)f 



i=l 



(13) 



If fc = —1 then E{x{t)) is constant along trajectories, 
however it is strictly decreasing for fc > —1 provided 
x{t) is not an equilibrium of (4) (i.e. whenever at least 
one of the phase differences a;, does not satisfy r(a;i) = 
r((5)). With this functional E and a bound on its rate 
of decrease (13) wc establish in Appendix A that the 
linearization J = Cr'(a;*) is hyperbolic for fc > — 1. This 
implies that eigenvalues of J cannot cross the imaginary 
axis when fc is varied. Thus, for all fc > — 1 the indices 
i>u{x^) and Vsix*) are the same as for fc = 0. Since for 
fc = the Jacobian J is an upper diagonal matrix with 
its eigenvalues on its diagonal, the number of negative 
diagonal entries of J is fo(a;,) (the diagonal entry of C 
is —1), and the number of positive diagonal entries of J 
is f7r(a;*). □ 
We observe that for (5 = the functional E{x) is 
bounded, leading for fc 7^ 1 to the result that full synchro- 
nization (the equilibrium (0, . . . , 0) for (5 = 0) is globally 
stable in the same way as the classical results [16]. 



V. ROTATING WAVES AND SADDLE 
CONNECTIONS 

Despite the existence of a Lyapunov functional the 
dynamics of (4) is not necessarily trivial because the 
phase space is an TV-dimensional torus. Gencrically we 
can expect that the unstable manifold of an equilib- 
rium x and the stable manifold of an equilibrium y in- 
tersect, giving rise to heteroclinic saddle connections, if 
i^TTix) + i'o{y) > N. This implies that for any equilibrium 
X heteroclinic connections to all equilibria y satisfying 
Ujriy) < i^Trix) are generic. 

Saddle connections between equilibria of the same 
type, that is ly-jrix) = v^{y), are of co-dimension 1 and 
can be achieved by tuning the parameters 5 and fc. These 
non-generic saddle connections also form co-dimension 1 
boundaries of periodic orbits of rotating wave type. 
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FIG. 2: Co-dimension 1 homoclinic connections and rotating 
waves for N = 2, r{x) = sin a: and 5 = 1. (a): bifurcation di- 
agram for rotating wave solutions (computed with (r)AUTO 
[17, 18]), (b)-(d): phase portraits showing equilibria homo- 
clinic connections and rotating waves for parameter values of 
k indicated in panel (a). The direction of the flow is always 
upward. 
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FIG. 3: Bifurcation diagram in the (fc, 5)-plane for N = 2 and 
V{x) = sin I (computed with (r)AUTO [17, 18]). 



Figure 2(a) shows for = 2 that periodic orbits (of 
rotating wave type) are possible for a range of param- 
eters k and 5 — \ despite the existence of a Lyapunov 
functional E. In the parameter space the existence of 
periodic orbits is bounded by two co-dimension 1 bifur- 
cations: a saddle-node of periodic orbits corresponding 
to the phase portrait in Figure 2(b) at the minimal k 
supporting periodic orbits and a heteroclinic connection 
between neighboring saddles (the a;2-component has in- 
creased by 27r along the heteroclinic connection). 

If the phase shift by 27r is ignored then the saddle con- 
nection is homoclinic. Figure 3 shows the parameter re- 
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FIG. 4: Bifurcation diagrams for A'^ = 3 and r(a;) = sin a:: 
(computed with (r)AUTO [17, 18]). Panel (a): family of 
periodic orbits for fixed fe = —0.5 and varying & and its 
codimension-1 bifurcations. Panel (b): region with rotating 
waves in the (fc, 5)-plane with its codimension-1 bifurcations. 



gion in the (A:, ^)-plane where periodic orbits of rotating 
wave type exist. The special points in the plane are at 
5 — 7^/2 when all four equilibria collapse into a single 
degenerate equilibrium and at fc = — 1 when the sys- 
tem is a single-degree-of-freedom conservative oscillator 
for TV = 2. We also note that the diagram is reflection 
symmetric about the line 5 = n /2 as expected. 

In dimensions N > 2 the dynamics can become chaotic. 
Numerical evidence for this is shown in Figure 4 for 
A'' = 3. Panel (a) shows a family of periodic orbits of ro- 
tating wave type {xi{T) , X2{T) , xsiT)) = (a;i(0), 2:2(0) + 
27r, 2:3(0) -I- 2tt) where T is the period, 5 is varied, and 
k — —0.5 is kept fixed. We observe that this family un- 
dergoes a period doubling sequence. Moreover, the lower 
(predominantly unstable) part of the branch reaches a 
homoclinic connection to the saddle Xg ^ {it — S, S, 5) 
at Hg. The eigenvalues at this saddle have the form 
{fi+^—^j.^ ± iu;_) where fj,± > and a;_ > are real 
numbers and /i_ < /i+ < 2/i_. Shil'nikov's results im- 
ply that there is an infinite number of period doubling 
cascades of stable periodic orbits close to the homoclinic 
connection Hg (under certain non-degeneracy conditions, 
see [19, 20]). The precise sequence of period-n branches 



5 



60 
40 

^3 

20 




60 









- 1 








/ 5 X 



10 



FIG. 5: Rotating periodic orbit of period 8 for = 3, r(a;) = 
sin a;, k = —0.5, 5 ~ 1.19, near the end of the period doubling 
sequence in parameter space. 



for a homoclinic to a saddle of this type, and how to cal- 
culate them numerically can be found in [21]. The other 
end of the family of periodic orbits is also a homoclinic 
connection, to the saddle ((5, tt — (5, (5), which has three real 
eigenvalues > > — > — /i-,2 where /i+ > 
This implies that there is only a single periodic orbit close 
to this homoclinic connection (no snaking), and that this 
periodic orbit is unstable [20]. 

In the two-parameter plane (shown in Figure 4(b)) we 
observe a shape that looks superficially similar to the case 
N = 2 except that the range of S over which the region 
of rotating waves extend is smaller. However, as one can 
see from Figure 4(a) the rotating waves are not stable 
inside the region bounded by the homoclinic connection 
Hg and the saddle-node. The period doubling sequence 
(also shown in Figure 4(a)) gives a better estimate for the 
region of stable periodic rotating waves. The remainder 
of the region is not filled with chaotic rotations because 
the chaotic attractor at the end of the primary period 
doubling cascade collides with heteroclinic saddle con- 
nections. Figure 5 shows a periodic motion of period 8, 
which is close to the end of the period doubling sequence 
in parameter space. 

Special points in the bifurcation diagram 4(b) are on 
the symmetry line 6 — 7^/2 (where all equilibria collapse 
to a single degenerate equilibrium) and various degen- 
eracies of the homoclinic connection, for example, the 
heteroclinic connection (tt — 6, 6, S) — {2tt + 6,tt — 6,6) — 
(tt — 5, 2Tr + 6, 2tt + 6) indicated in Figure 4(b) (this degen- 
eracy ends the numerical continuation of the homoclinic) . 

In the original variables 9j the rotating waves shown 
in this Section correspond to a continually drifting phase 
difference between the oscillators 2 and 3 (for N — 2) 
while the phase difference between oscillators 1 and 2 re- 
mains bounded. The homoclinic boundary corresponds 
to regimes where two neighboring oscillators hover in 
anti-phase for a long time before a phase slip occurs. 



For A'' = 3 the chaotic regimes near the Shil'nikov saddle 
correspond to the regime where oscillators 1 and 2 are 
nearly in anti-phase and 2, 3 and 4 are nearly in-phase 
for a long time before the phase differences between 2 
and 3, and 3 and 4 slip (nearly) simultaneously. 
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FIG. 6: Fraction of phase and parameter space not attracted 
to synchronization. These are simulation results obtained 
during 100 trials where in each trial we chose 5, 000 ran- 
dom initial conditions a;(0) € S'^ and random parameters 
S e [0,7r/2], k e [-1,0] for small N > 2. During each trial 
simulation we record when a trajectory enters the neighbor- 
hood of the stable equilibrium {S, . . . ,S). We stop the sim- 
ulation if no new trajectory enters this neighborhood during 
a time period of 200. The plot (horizontal line) shows which 
fraction of the 5, 000 realizations was not yet trapped on av- 
erage when we stopped the simulation. The vertical lines give 
the variability between trials (25% and 75% quartile) and the 
crosses are the extreme (maximal and minimal) result of the 
100 trials. 

As the number of oscillators is increased by 1, we ob- 
serve that the parameter region which supports stable 
periodic rotations shrinks. Put another way, full syn- 
chronization becomes more prevalent for iV = 3. The 
simulation results shown in Figure 6 support this obser- 
vation. This is in contrast to the common feature that 
longer chains of oscillators require stronger coupling for 
synchronization [16, 22] in symmetrically coupled chains 
with random frequencies. The reason behind this ap- 
parent difference is that in our setup the only detuned 
oscillators are at the boundary (indices and N), and 
that the effect of the boundary diminishes for increasing 
N. We did not observe, however, that the fraction of 
slipping realizations approaches zero for large TV. 



VI. CONCLUSION 

In this paper we studied the dynamical system describ- 
ing the phase differences of a uniform chain of oscillators. 
Two important parameters (other than the number of el- 
ements in the chain) are the detuning of the end points, 
6 and the ratio k of the coupling strengths upward and 
downward the chain. We have classified the phase locked 
states by relating the number of phase differences that 
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equal tt — 5 to the number of unstable directions that 
this phase locked state has as a saddle in the phase space. 
We also studied periodic phase slipping (rotating waves) 
systematically for = 2 (3 oscillators) and A'' = 3 (4 os- 
cillators). The most curious difference between the case 
of N = 2 and N = 3 is the shrinking of the parameter re- 
gion where periodic rotations occur. Simulations suggest 
that their basin of attraction also shrinks as a fraction 
of the whole phase space volume. This would suggest 
that longer chains have more robust phase locking. This 
effect appears to hold with respect to both aspects of ro- 
bustness combined: the fraction of initial values (treating 
parameters also as dynamic variables with trivial dynam- 



ics, (5 = 0, fc = 0) in the phase space leading to locking is 
smallest for N = 2. 
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Appendix A: Proof of hyperbolicity of equilibria 

Proposition 2 (Absence of eigenvalues with zero real 
part). Let k > —1, and let x* be an equilibrium of (4). 
Then the Jacobian J — Cr'(x*) of (4) in cannot have 
eigenvalues on the imaginary axis. 
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Proof First wc note that the matrix J is regular for 
k > —1 (and, thus, cannot have an eigenvalue 0): 



detC 



_ (-1)^ (fc^+i - 1) 

detr'(x*) = r'((5)^(-i)"' 



7^0, 

^0. 



(Al) 



We show the absence of purely imaginary eigenvalues in- 
directly. In preparation for this part of the proof we es- 
tablish a bound on the gradient of Lyapunov fimctional 
E, given in (12). In an equilibrium a;* the gradient of 
E vanishes and the Hessian of E in a;* is equal to 
r'(a;.): 
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;E{X%=X, 



{dxY' 
= p ■ diag((Ti, . . . ,o-„ 



(A2) 



r'(x*). 



Consequently, the equilibrium (5,. . . ,5) is a local mini- 
mum of E, and (tt— (5, . . . , it — 5) is a local maximum of E. 
Equilibria which have some components equal to 5 and 
others equal to it — 5 are saddle points of the graph of E. 
As we want to study the local stability of an equilibrium 
a;* we introduce quantities measuring the deviation from 

y = x - X* and D{y) = E{y + a;*) - ^^(a:*), (A3) 

such that D{0) = and 

Diy)=y'^H.y + 0{\\yf), (A4) 

for all 2/ in a neighborhood of 0. Furthermore, equation 
(13) estimates D from above for k > —1: 



where cq > is a constant independent of y. 

Assume that J has a purely imaginary eigenvalue, say, 
X = ijjb where /x > 0, and let u be the corresponding 
eigenvector. We choose the time T = 2n/fi. Let c > 
be sufficiently small (e will depend on T). After time T 
the solution of (4) starting from a;(0) = a;* + eu (let us 
call the solution a;(-)) satisfies x{T) = a;, + eu + 0{e'^). 
Consequently, 

D{x{T)) - D{x{0)) = 

= [a:(T) - x,\^H,\x{T) ~ x.] + 0(||a;(T) - x^^) 

- e^u^H^u - 0{e^) 
= O(e^). (A6) 

On the other hand, the trajectory x{t) for t G [0,T] is 
a perturbation of order of an ellipse with a minimal 
radius of order e around x,. This means that we can 
choose a uniform constant ci of order 1 such that cie is 
smaller than this minimal radius (ci is uniform in e) and, 
hence, 

\\x(t)~x4>cie (A7) 
for all t e [0,T]. Inequality (A5) implies that 

D{x{T)) - D{x{Q)) = f b{x{t))dt 
Jo 

—Cq I \\x{t) — Xt:\f dt (due to (A5)), 

(due to (A7)), 



< 



,2,2 



D{y) < -co\\y\\ 



(A5) 



< -coTcfe 



which contradicts (A6) because the constants cq and ci 
are independent of e. □ 



